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ABSTRACT 

A  modeling  algorithm  is  presented  to  compute  simultaneously  polymer  con¬ 
formations  and  ionic  current,  as  single  polymer  molecules  undergo  translo¬ 
cation  through  protein  channels.  The  method  is  based  on  a  combination  of 
Langevin  dynamics  for  coarse-grained  models  of  polymers  and  the  Poisson- 
Nernst-Planck  formalism  for  ionic  current.  For  the  illustrative  example  of 
ssDNA  passing  through  the  a-hemolysin,  vivid  details  of  conformational 
fluctuations  of  the  polymer  inside  the  vestibule  and  /5-barrel  compartments 
of  the  protein  pore,  and  their  consequent  effects  on  the  translocation  time 
and  extent  of  blocked  ionic  current  are  presented.  In  addition  to  shedding 
insights  into  several  experimentally  reported  puzzles,  our  simulations  offer 
experimental  strategies  to  sequence  polymers  more  efficiently. 
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1  Introduction 


Translocation  of  polymers  through  biological  channels  is  very  complex  in¬ 
volving  many  machineries,  and  is  a  fundamental  step  in  many  life  processes. 
While  several  essential  features  of  translocation  are  richly  documented  in 
systems  such  as  mRNP  complex  through  nuclear  pores  [1-3],  only  recently 
a  simple  system  has  been  identified  to  follow  the  single-file  passage  of  one 
isolated  polymer  through  one  channel  [4-11].  In  this  system,  the  channel 
is  constituted  by  self- assembling  heptamers  of  the  Staphylococcus  aureus 
cr-hemolysin  (aHL)  protein.  The  channel  is  assembled  in  a  phospholipid 
bilayer,  which  offers  a  physical  barrier,  and  the  channel  has  an  opening  of  di¬ 
ameter  of  about  1.4  nm  at  the  narrowest  constriction  [12].  A  single-stranded 
polynucleotide,  such  as  poly(deoxyadenylate)  and  poly(deoxycytidylate),  is 
pulled  through  the  channel  by  an  externally  applied  voltage  gradient  across 
the  channel  in  a  solution  of  a  strong  electrolyte.  The  idea  is  that  the  ionic 
current  through  the  channel  due  to  the  passage  of  small  ions  of  the  electrolyte 
is  blocked  to  certain  extent  during  the  event  of  translocation  of  the  polymer. 
It  has  been  hoped  that  the  extent  and  duration  of  the  current  blockade  are 
unique  signatures  of  the  identity  of  the  polymer,  both  in  terms  of  polymer’s 
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chemical  characteristics  and  physical  length. 

Even  this  simplest  setup,  where  identical  molecules  undergo  transloca¬ 
tion,  has  generated  several  puzzling  results.  The  distribution,  P(r),  of  the 
duration  r  of  blockade  of  ionic  current  lb  is  very  broad  and  appears  to  exhibit 
at  least  two  peaks.  In  addition,  there  are  also  several  levels  of  ionic  current 
blockade  for  the  same  molecule.  It  is  a  standard  practice  in  experimental 
investigations  to  combine  the  histograms  of  r  and  7ft  [8] .  The  resultant  scatter 
plots  always  yield  two  groups  of  data  even  for  monodisperse  homopolymers. 

To  gain  insight  into  these  puzzles,  we  have  developed  the  following  simu¬ 
lation.  This  is  complementary  to  a  flurry  of  theoretical  activity [13-21],  based 
on  entropic  barrier  dynamics[22],  all  of  which  leading  to  a  generic  P(r)  unlike 
in  experiments.  Although  it  is  indeed  desirable  to  perform  the  computation 
ab  initio ,  the  size  of  the  system  to  be  simulated  is  forbiddingly  large  to  enable 
such  a  computation.  Therefore  we  restrict  ourselves  to  generate  a  minimal 
model  of  polymer  translocation  through  aHL  channel,  and  our  goal  is  to 
discover  the  most  basic  concepts  relevant  to  a  molecular  understanding  of 
the  experimental  results.  Our  minimal  model  incorporates  enough  chemi¬ 
cal  details  of  the  polymer  and  channel  to  evaluate  the  potential  roles  played 
by  secondary  structures  of  the  polymer,  orientation  of  the  chain  (3’  versus 
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5’  end),  binding  sites  inside  the  lumen  of  the  channel,  etc.  In  addition  to 
resolving  the  above  puzzles,  the  details  emerging  from  our  simulations  are 
so  vivid  that  new  experimental  protocols  can  be  formulated  to  substantially 
narrow  down  P(t),  and  thus  to  enable  a  much  faster  sequencing  scheme. 

The  key  result  of  our  simulations  is  that  conformational  entropy  of  the 
polymer  plays  a  dominant  role  in  dictating  r  and  P(r).  The  breadth  of  P(r) 
can  be  directly  attributed  to  the  entropic  trap  arising  from  the  vestibule  of 
oHL  pore.  P(r)  can  be  made  dramatically  sharper,  sufficient  to  facilitate 
sequencing  of  polymers,  by  deleting  the  vestibule  or  by  dragging  the  polymer 
in  a  direction  normal  to  the  pore  axis. 

2  Simulation  Method 

The  system  consists  of  essentially  four  parts,  namely,  the  heptameric  a- 
hemolysin  pore,  a  ssDNA  molecule  undergoing  translocation,  a  phospholipid 
membrane  barrier  carrying  the  pore,  and  an  electrolyte  solution  permeating 
through  the  pore.  The  protein  pore,  aHL,  is  represented  by  a  united  atom 
model  where  the  residues  of  the  pore  are  replaced  by  beads.  The  coordinates 
of  the  residues  are  extracted  from  the  Protein  Data  Bank.  There  can  be  many 
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choices  for  the  description  of  beads  in  terms  of  their  sizes  and  shapes.  In  this 
paper  we  present  results  for  the  simplest  choice  of  the  beads.  In  the  present 
united  atom  model,  each  residue  is  replaced  by  a  spherical  bead  and  the 
diameter  of  each  bead  is  3 A.  A  bead  is  assigned  -1,  0,  or  +1  charge  according 
to  the  ionization  of  the  residue.  All  aspartic  acid  and  glutamic  acid  residues, 
and  C-terminal  groups  are  deprotonated;  arginine,  lysine,  histidine  and  N- 
terminal  groups  are  protonated.  The  rest  of  the  residues  are  electrically 
neutral.  In  view  of  experimental  contexts  [23],  the  protein  pore  is  not  allowed 
to  move. 

Further,  the  protein  pore  is  treated  as  a  rigid  structure  with  a  dielectric 
constant  of  2  surrounded  by  a  solution  of  dielectric  constant  80.  The  protein 
pore  is  embedded  in  a  48 A  thick  membrane  of  dielectric  constant  2.  The 
membrane  surfaces  are  hard  walls.  The  protein  pore  is  the  only  space  to 
allow  passage  of  the  polymer  from  the  donor  compartment  to  the  recipient 
compartment.  The  symmetric  axis  of  the  pore  is  positioned  along  the  X-axis, 
as  shown  in  Fig.  (la)  which  gives  the  geometry  of  the  pore.  The  vestibule 
of  the  pore  is  submerged  in  the  donor  compartment.  The  /5-barrel  of  the 
protein  begins  at  x  =  0  and  extends  up  to  x  =  48A. 

We  have  used  the  side-chain  incorporated  model  to  represent  a  single- 
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stranded  polynucleotide  [20],  due  to  its  simplicity,  fast  computational  speed 
and  past  success.  Each  nucleotide  is  made  up  of  three  subunits  (phosphate, 
sugar,  and  base)  and  each  subunit  is  taken  as  a  bead.  The  bead  representing 
the  phosphate  moiety  carries  one  negative  charge.  To  distinguish  the  chain 
directionality  in  terms  of  3’  versus  5’  ends,  we  deliberately  tilt  the  side-chain 
(base)  at  65°  angle  as  shown  in  Fig.  (lc).  The  chain  backbone  is  freely  jointed 
in  view  of  the  fact  that  the  persistence  length  of  single-stranded  DNA  has 
been  measured  in  high  salt  concentrations  to  be  roughly  the  length  of  1-2 
nucleotides  [24-26].  Although  the  three  subunits  are  of  different  sizes  and 
masses,  the  results  presented  below  are  for  the  model  of  all  subunits  having 
equal  size  and  mass.  The  diameter  of  these  beads  is  taken  to  be  2.5A,  and 
the  equilibrium  “bond  length”  between  these  “atoms”  is  2. 5 A.  The  mass  of 
a  nucleotide  is  distributed  equally  into  the  three  subunits.  For  example,  for 
the  case  of  poly(dC),  each  subunit  carries  96Da. 

Polymer  conformations  are  evolved  by  using  the  velocity  Verlet  algorithm 
[27]  to  follow  the  dynamics  of  the  ffh  bead, 

m^  =  -Cvi-VriU  +  Fi(t)  (2.1) 
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where  r),  m,  and  C  are  the  position  vector,  mass,  and  friction  coefficient, 
respectively,  of  the  ith  bead.  F(t)  is  the  random  force  from  the  solvent  bath 
acting  on  the  ith  bead  and  is  stipulated  to  satisfy  fluctuation-dissipation 
theorem, 

<  Fi{t)  ■  F3(t')  >=  5i:j6kBT(6(t  -  t')  (2.2) 

where  t  is  the  time  and  5  is  the  usual  delta  function.  U  in  Eq.(2.1)  is  the  total 
potential  energy  acting  on  the  ith  bead,  and  consists  of  four  contributions. 

U  =  Ubond  +  Ulj  +  Udh  +  Uy-  (2-3) 

These  correspond,  respectively,  to  bond  stretch,  short  range,  screened  Couloumb, 


and  electric  potential  given  by 

^bond  =  ~  M2  (2.4) 

Ulj  =  eLj[(°r  ~  2(^)6]  (2.5) 

z  z 

UDH(rij)  =  4^  ^rexP(~^r)  (2.6) 

Uv{?i)  =  ZiVifi).  (2.7) 
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Here,  l  is  the  bond  length,  l0{=  2.hA)  is  the  equilibrium  bond  length,  r  is  the 
distance  between  interacting  beads,  Zi  is  the  electric  charge  on  the  ith  bead, 
eo  is  the  permittivity  of  the  vacuum,  and  V  is  electric  potential.  The  values 
of  the  parameters  are:  k  =  171  heal /mol  ■  A2,  tLJ  =  0.5  heal /mol,  cr  =  2.5A, 
k~1  =  3\  for  1  M  KC1.  The  dielectric  constant  e  is  inhomogeneous,  e  is  2 
inside  beads  and  membrane  and  is  80  everywhere  else. 

The  ionic  current  I(t )  at  time  £,  due  to  the  passage  of  electrolyte  ions 
accompanying  the  polymer  transport,  is  computed  by  using  the  Poisson- 
Nernst-Planck  (PNP)  formalism  [28-30].  Taking  advantage  of  the  fact  that 
small  ions  relax  much  faster  than  a  large  polyelectrolyte  molecule  and  that 
the  concentration  of  electrolyte  in  the  pertinent  experiments  is  very  high  in 
comparison  with  monomer  concentration,  we  assume  that  at  every  time  step 
of  the  Langevin  dynamics  simulation  of  the  polymer,  the  electrolyte  ions  have 
relaxed  to  the  steady  state  so  that  the  polymer  chain  is  taken  only  as  a  fixed 
charge  distribution  pp(r,  f)  at  this  time. 

With  this  assumption,  the  self-conSistent  coupled  PNP  equations  for  the 
local  charge  density  C,(r,  t)  of  the  ith  ionic  species  and  the  electric  potential 
V  (r,  t)  for  a  given  polymer  conformation  at  t  are 
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v '  (vci(r)  +  f~f  ww) =  °- 


(2.8) 


e0V-[c(r)W(r)]  -  -p/(r)  -  pp(r,t)  -  ^Ci(r).  (2.9) 

i 

Here  p/(r)  is  the  local  charge  density  arising  from  fixed  charges  in  the 
system,  such  as  the  charged  beads  of  the  protein  pore.  As  already  mentioned, 
e(r)  is  set  to  2  for  the  beads  of  polymer  and  pore,  and  for  the  membrane.  It 
is  80  everywhere  else.  In  addition,  the  externally  applied  voltage  gradient  is 
accounted  for  by  fixing  the  electrostatic  potential  at  the  boundaries  in  the 
x-direction.  The  box  for  the  ionic  current  calculations  is  located  at  -100  A  < 
X  <  75  A,  -50  A  <  Y  <  50  A,  and  -50  A  <  Z  <  50  A.  We  set  V  —  0  at 
X  =  -100A  and  V  =  Fo  at  A  =  75A. 

The  above  equations  (2.8)  and  (2.9)  are  solved  by  successive  over-relaxation 
method  with  a  grid  spacing  of  lA,  by  following  Ref. (29).  First  we  start  with 
V(X)  due  to  the  externally  applied  electrostatic  potential  gradient  across  the 
membrane  of  dielectric  constant  of  2,  and  uniform  concentrations  (C+  and 
CL)  of  cations  and  anions  of  the  dissolved  electrolyte. 

Then,  the  coordinates  of  beads  of  the  pore  and  polymer  are  read  in.  Next, 
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using  the  successive  over-relaxation  method,  V,  C+,  and  CL  are  computed 
iteratively  until  they  converge  within  the  tolerance  levels  of  AF  —  10  8 V, 
and  A C±  =  10 ~7M.  In  our  simulations  we  have  considered  1  M  of  KC1  as 
the  electrolyte.  Since  C±  is  very  high  in  comparison  with  the  counterions 
from  the  polymer,  the  latter  are  ignored.  From  the  convergent  values  of  C± 
and  V  at  each  time  unit  of  the  Langevin  simulation,  the  ionic  current  at  this 
time  t  is  given  by 

I(t)  =  A(J+  +  J-),  (2.10) 

with 

J±M  =  -D±  ( VC±(r,t )  +  ^W(r,t))  .  (2.11) 

In  evaluating  the  ionic  current,  we  have  taken  J±  at  X  —  48A  and  A  at 
this  exit  point  is  44lA2.  D+  and  £>_  are  taken  to  be  1.96  x  10~5  cm2s_1  and 
2.03  x  10“5  cm2s_1,  respectively,  for  K+  and  Cl-  ions.  The  magnitude  of  the 
applied  voltage  is  Vo  =  120mVL 

The  simulation  is  carried  out  as  follows.  First,  a  single  polymer  chain  is 
equilibrated  in  the  absence  of  pore  and  externally  applied  voltage  gradient,  by 
considering  only  the  connectivity,  Lennard-Jones,  and  Debye-Hiickel  forces 
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among  various  beads  of  the  chain.  Separately  the  electrostatic  potential  for 
the  protein  pore  is  computed  for  1  M  KC1  solution  in  the  steady  state  with  a 
given  value  of  Vo (120 mV),  and  in  the  absence  of  the  polymer.  The  potential 
along  the  central  axis  of  the  channel  is  given  as  dotted  numerical  data  in  Fig. 
(2) .  The  different  gradients  (slopes  of  the  curve)  at  different  parts  of  the  tra¬ 
jectory  are  readily  conspicuous  in  this  figure.  Since  we  are  interested  here  in 
the  details  of  the  actual  translocation  event,  and  not  in  the  way  the  polymer 
arrives  at  the  channel  entrance,  we  have  extended  the  range  of  the  potential 
gradient  inside  the  vestibule  to  the  donor  compartment.  In  addition,  we  have 
ignored  the  finer  details  of  the  potential  gradient  and  replaced  the  actual  po¬ 
tential  gradient  by  the  approximated  gradient  consisting  of  three  parts,  viz., 
up  to  the  /3-barrel,  inside  the  /3-barrel,  and  outside  /3-barrel  in  the  recipient 
compartment. 

Next,  the  equilibrated  polymer  is  inserted  in  the  donor  compartment  such 
that  the  center  of  mass  of  the  chain  is  at  X  =  —  100A  in  front  of  the  mouth 
of  the  pore.  The  polymer  is  dragged  by  the  electrical  potential  gradient 
across  the  pore.  For  each  time  unit  of  the  Langevin  dynamics  simulation, 
the  coordinates  of  the  polymer  beads  are  stored  and  the  ionic  current  at 
this  time  is  computed.  When  any  one  monomer  of  the  polymer  enters  the 
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mouth  of  the  pore  (at  X  =  —50 A)  a  clock  for  measuring  the  translocation 
time  starts.  As  the  polymer  progressively  invades  the  lumen  of  the  pore,  the 
ionic  current  is  calculated  at  each  time  unit.  This  is  how  the  ionic  current 
traces  given  below  are  constructed.  It  is  to  be  noted  that  VV  entering  the 
Langevin  equation  for  the  polymer  bead  is  the  original  value  of  the  protein 
pore  in  the  absence  of  the  polymer.  When  the  last  bead  of  the  polymer  exits 
the  outer  edge  of  the  pore  (at  X  =  48 A),  the  clock  stops  with  a  reading  of 
the  translocation  time  r.  By  repeating  the  above  sequence  of  simulations 
thousands  of  times  we  construct  the  histogram  of  r,  and  for  each  of  the 
entries  of  the  histogram,  we  have  a  trace  of  ionic  current  as  the  polymer 
traverses  through  the  pore. 

3  Results 

First,  we  give  results  on  chain  conformations,  ionic  current,  and  translocation 
time  t  as  the  polymer  transits  through  the  channel  in  typical  simulations. 
Then  we  present  distribution  functions  of  blocked  current  and  r,  based  on 
thousands  of  simulations.  Fig.  (3a)  represents  a  typical  simulation  result 
for  the  model  chain  of  ss-poly(dC)  with  N  =  45  nucleotides  undergoing 
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translocation  through  crHL  under  a  potential  difference  of  120  mV  in  1  M 
KC1  solution.  One  unit  of  simulation  time  is  arbitrarily  taken  to  be  0.15  [is, 
by  matching  the  peak  values  of  r  in  our  simulations  and  the  experiments  of 
Refs.  (8)  and  (9),  and  by  noting  that  the  average  translocation  velocity  is 
constant  in  the  long  polymer  and  high  voltage  gradient  limits  [9]. 

There  are  several  key  features  evident  from  the  typical  trajectory  of  Fig. 
(3a).  (i)  The  open  pore  current  of  about  135  pA  is  very  close  to  the  ex¬ 
perimental  value,  (ii)  The  polymer  enters  the  mouth  of  the  vestibule  not 
necessarily  as  a  single-file.  The  presence  of  polymer  monomers  at  the  mouth 
of  the  vestibule  is  sufficient  to  provoke  a  response  in  the  ionic  current  by 
providing  a  spatial  blockade  to  the  flow  of  electrolyte  ions,  (iii)  As  time 
progresses  the  chain  gets  into  the  vestibule.  Due  to  its  large  volume,  the 
vestibule  acts  as  an  entropic  trap  and  the  chain  segments  linger  in  this  cav¬ 
ity  for  some  definite  time  duration.  At  this  stage,  the  ionic  current  is  If, 1  (~ 
100  pA).  The  fluctuations  in  I  around  its  mean  value  reflect  the  dynamics 
of  the  polymer  rattling  inside  the  vestibule,  (iv)  One  end  of  the  chain  even¬ 
tually  enters  the  /3-barrel.  At  this  time,  I(t )  is  reduced  precipitously  to  the 
blocked  current  level  Ib 2.  The  average  value  of  Ib2  is  17  pA.  This  very  low 
value  reflects  the  reduction  in  the  cross-sectional  area  of  /3-barrel  for  flow  of 
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electrolyte  ions  by  the  presence  of  polymer  segments.  The  fluctuations  in  h'2 
is  solely  due  to  the  dynamics  of  the  polymer.  The  ionic  current  is  about  J52 
until  the  last  monomer  exits  the  end  of  the  /5-barrel,  r  is  the  duration  of 
all  of  the  above  events,  (v)  When  the  chain  has  completed  its  translocation, 
the  ionic  current  returns  to  the  open  pore  current  after  leaving  a  distinct 
signature  for  the  presence  of  some  monomers  in  the  back  of  the  exit  location, 
perturbing  the  flow  of  electrolyte  ions. 

It  is  to  be  noted  from  Fig.  (3a)  that  the  chain  spends  significant  amount  of 
time  inside  the  vestibule  before  entering  the  /5-barrel.  While  such  trajectories 
are  quite  common,  there  are  also  more  prevalent  trajectories  where  the  dwell 
time  inside  the  vestibule  is  short  as  illustrated  in  Fig.  (3b).  The  essential 
features  of  Fig.  (3b)  are  the  same  five  features  discussed  for  Fig.  (3a), 
except  that  the  time  spent  inside  the  vestibule  is  shorter.  As  seen  from  the 
chain  conformations  presented  in  Fig.  (3b),  the  shorter  duration  inside  the 
vestibule  is  due  to  the  initial  chain  orientation  in  line  with  the  axis  of  the 
channel,  and  the  chain  entering  the  vestibule  mouth  essentially  as  a  single¬ 
file. 

Vivid  details  such  as  those  in  Figs.  (3a)  and  (3b)  are  available  for  each  of 
the  thousands  of  simulations  that  were  performed,  based  on  which  we  have 


14 


constructed  histograms  for  r  and  /*,.  The  normalized  probability  P(r)  based 
on  4000  simulations  is  plotted  against  r  in  Fig.  (4)  for  N  =  45,  Vo  =  120 mV, 
and  1  M  KC1.  The  error  bars  are  estimated  from  the  square  root  of  counts  per 
bin  before  normalizing  the  histogram.  It  is  evident  that  the  distribution  of  r 
is  very  broad  and  two  dominant  peaks  at  50  (is  and  75  (is  may  be  identified. 
This  result  is  in  agreement  with  our  earlier  report  of  two  peaks  based  on 
smaller  number  of  simulations  with  shorter  (N  =  15)  chains  [20].  Now  we 
are  in  a  position  to  go  back  to  the  trajectories  of  the  chains  and  study  them 
for  various  events  taking  different  translocation  times.  Indeed,  an  example 
of  the  trajectory  taking  76  (is  is  Fig.  (3a),  where  the  vestibule  plays  the 
role  of  an  entropic  trap,  delaying  the  translocation  process.  In  contrast,  by 
stochasticity  associated  with  chain  end  orientation  and  thermal  noise  from 
the  background,  the  chain  end  can  enter  the  /3-barrel  sooner  by  essentially 
avoiding  the  vestibule’s  entropic  trap.  This  is  the  case  with  Fig.  (3b)  where 
t(=  23.7 (is)  is  much  shorter,  and  this  event  belongs  to  the  dominant  peak  in 
the  histogram.  We  validate  this  mechanism  for  the  occurrence  of  two  peaks, 
by  performing  further  simulations  on  nanotubes  where  we  have  cut  out  the 
vestibule  part  of  the  channel  (see  below).  We  must  mention  at  this  point 
that  there  is  no  significant  discrimination  between  3’  end  and  5’  end  entering 
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the  ^-barrel  first  in  terms  of  translocation  time. 

As  already  noted,  there  are  two  levels  of  blocked  ionic  current  during  the 
translocation  event.  For  the  event  corresponding  to  Fig.  (3a),  the  normal¬ 
ized  histogram  P(h)  of  blocked  current  It,  within  the  duration  of  r,  given  in 
Fig.  (5a),  exhibits  two  separate  distributions.  The  average  values  of  these 
two  peaks  represent  whether  the  blockade  corresponds  to  the  vestibule  (weak 
blockade)  or  the  /3-barrel  (strong  blockade).  The  widths  of  the  peaks  corre¬ 
spond  to  the  accompanying  polymer  dynamics.  It  must  be  remarked  that  the 
occurrence  of  two  populations  of  blocked  ionic  current  has  nothing  to  do  with 
the  actual  value  of  r.  The  histogram  of  Fig.  (5a)  corresponds  to  r  =  76 [is. 
For  the  case  of  r  =  23.7 fxs  corresponding  to  Fig.  (3b),  the  histogram  of  Ib 
is  given  in  Fig.  (5b).  Here  again,  there  are  two  distributions  representing 
polymer  dynamics  inside  the  vestibule  and  the  /3-barrel.  However  the  weight 
of  A  population  relating  to  the  vestibule  is  weaker  for  the  faster  transloca¬ 
tion  events.  These  results  emphasize  the  need  to  interpret  the  scatter  plots  of 
Ref.  [8]  differently  by  having  to  go  into  individual  events  instead  of  clumping 
them  all  together. 

The  intervention  by  the  vestibule  of  aHL  as  an  entropic  trap  for  the 
translocation  of  the  polymer  leads  to  the  breadth  of  P(r).  To  further  validate 
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this  observation,  we  have  repeated  our  simulations  by  replacing  the  aHL 
channel  by  a  cylindrical  tube  of  diameter  18  A  and  length  98  A,  as  illustrated 
in  Fig.  (lb).  The  tube  is  inserted  inside  the  membrane  of  dielectric  constant 
2,  and  the  membrane  thickness  is  the  same  as  the  tube  length. 

A  typical  trajectory  of  polymer  translocation  through  the  cylindrical  tube 
is  given  in  Fig.  (6)  for  N  =  45,  Vo  =  120 mV,  and  1  M  KC1.  This  also 
exhibits  all  of  the  features  described  for  Figs.  (3a)  and  (3b),  except  the 
part  corresponding  to  the  role  of  the  vestibule,  where  the  weaker  blockade 
is  prominently  absent.  Based  on  1000  simulations,  the  histogram  P(t)  is 
given  in  Fig.  (7),  where  the  results  of  Fig.  (4)  are  included  for  comparison. 
It  is  clear  that  P(r)  is  very  sharply  peaked  and  there  is  only  one  peak  for 
the  nanotube.  Our  simulations  thus  suggest  that  it  is  preferable  to  have  a 
nanotube  or  just  the  /9-barrel  for  sequencing  purposes  instead  of  aHL. 

4  Discussion  and  conclusions 

Our  simulations  using  coarse-grained  models  of  the  polymer  and  the  channel 
with  appropriate  accounting  of  dielectric  heterogeneity  reproduce  all  of  the 
essential  features  observed  experimentally  with  aHL.  The  vestibule  is  found 
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to  play  a  significant  role  by  providing  an  entropic  trap  and  thereby  slowing 
down  the  process.  Another  observation  is  that  there  is  no  direct  correlation 
between  translocation  time  and  blocked  current.  For  each  value  of  r,  there 
are  two  dominant  values  of  blocked  current. 

As  already  mentioned,  there  is  a  perturbation  in  the  ionic  current  through 
the  channel,  as  soon  as  a  monomer  is  eclipsing  the  pathway  for  small  ions  in 
front  of  the  channel  (see  Figs.  3  and  6).  This  suggests  a  novel  experimental 
setup  for  sequencing  polymers  using  measurements  of  ionic  current  through 
channels.  The  strategy  is  illustrated  in  Fig.  (8a),  where  the  circle  represents 
the  mouth  of  the  pore  through  which  electrolyte  ions  pass  and  create  the 
measured  ionic  current  under  an  applied  voltage  gradient.  The  pore  can 
be  either  ccHL  or  a  nanotube  with  prescribed  diameters.  Then  a  polymer  is 
dragged  across  the  front  of  the  pore  in  the  direction  of  the  arrow  in  Fig.  (8a). 
The  physical  size  of  the  monomer  at  the  pore-front  excludes  small  ions  to 
pass  by  and  consequently  reduces  the  magnitude  of  the  ionic  current.  When 
monomers  with  different  sizes  are  dragged  normal  to  the  ion  flow,  different 
levels  of  ionic  current  are  registered.  When  a  polymer  with  eight  dC  units 
is  dragged  at  the  speed  of  1.0 A//j,s  near  the  pore  entrance,  in  the  normal 
direction  to  the  current  flow,  the  time-dependence  of  ionic  current  is  given 
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in  Figs.  (8b)  and  (8c).  The  saw-tooth  nature  of  the  trace  is  due  to  the 
periodicity  arising  from  the  monomer  entering  the  aperture,  then  generating 
maximum  eclipse,  and  then  exiting  the  aperture.  If  an  octomer  of  dA  is 
simulated  under  identical  conditions  (where  the  diameter  of  the  side-chain 
bead  is  5 A  instead  of  2.5  for  dC  and  the  bond  length  connecting  the  side-chain 
to  the  backbone  is  3.75A),  the  ionic  trace  is  as  given  in  Figs.  (8b)  and  (8c).  In 
these  particular  simulations,  the  pore  diameter  and  length  are  taken  as  14A 
and  50A,  respectively,  to  detect  the  identity  of  only  one  monomer.  Depending 
on  the  size  of  the  aperture,  signals  corresponding  to  dimers,  trimers,  etc., 
can  be  generated.  In  the  present  illustration  of  the  concept,  we  now  drag 
octomers  with  different  sequences  containing  dA  and  dC.  Depending  on  which 
monomer  sequence  is  being  dragged,  the  ionic  current  will  toggle  back  and 
forth  between  the  values  corresponding  to  pure  dA  and  dC  values.  There  is 
only  one  unique  ionic  current  trace  for  a  given  sequence,  as  illustrated  by  two 
distinctly  different  traces  for  two  sequences.  Although  the  difference  in  the 
ionic  current  from  our  simulations  are  rather  small,  the  current  experimental 
status  [31]  is  able  to  discern  such  small  difference.  The  strategy  emerging 
from  our  simulations  will  enable  the  precise  sequencing  of  polymers  by  several 
orders  of  magnitude  faster  than  the  currently  available  techniques. 
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In  summary,  our  combination  of  Langevin  dynamics  simulations  of  poly¬ 
mer  dynamics  with  the  PNP  formalism  for  ionic  current,  in  the  general 
premise  of  united  atom  description,  offers  a  computational  platform  to  dis¬ 
cover  generic  physical  principles  behind  polymer  translocation  through  bio¬ 
logical  channels. 
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5  Figure  Captions 


Figure  1.  (a)  The  end  and  side  views  of  aHL.  (b)  The  internal  wall  of  the  nan¬ 
otube  is  made  of  spherical  neutral  beads  on  a  curved  hexagonal  lattice.  The 
interaction  between  polymer  bead  and  tube  bead  is  taken  as  the  Lennard- 
Jones  potential  of  Eq.  2.5  with  e  =  0.2  kcal/mol  and  a  =  2.75A.  (c)  United 
atom  representation  of  ss-poly(dC). 

Figure  2.  PNP  and  approximated  V(X)  across  the  protein  pore. 

Figure  3.  Typical  ionic  current  traces:  (a)  long  t,  (b)  short  r. 

Figure  4.  Histogram  of  r  for  aHL  (N=45,  V0=120  mV). 

Figure  5.  Histogram  of  /<,.  (a)  and  (b)  correspond,  respectively,  to  Figs.  3(a) 
and  3(b). 

Figure  6.  Typical  ionic  current  trace  for  a  nanotube. 

Figure  7.  Comparison  of  JP(t)  between  aHL  and  nanotube. 


25 


Figure  8.  (a)  Proposed  experimental  strategy.  The  ionic  current  traces  for 
(b)  two  homopolymers  and  (c)  two  heteropolymers  (sequence  given  in  the 
insert).  Each  polymer  has  8  bases. 
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